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£SJ , We present a first-principles derivation of the variational equations describing the dynamics of 

the interaction of a spatial soliton and a surface plasmon polariton (SPP) propagating along a 
metal/dielectric interface. The variational ansatz is based on the existence of solutions exhibiting 
differentiated and spatially resolvable localized soliton and SPP components. These states, referred 
to as soliplasmons, can be physically understood as bound states of a soliton and a SPP. Their 
respective dispersion relations permit the existence of a resonant interaction between them, as 
pointed out in The existence of soliplasmon states and their interesting nonlinear resonant 
behavior has been validated already by full-vector simulations of the nonlinear Maxwell's equations, 
' as reported in Q. Here, we provide the theoretical demonstration of the nonlinear resonator model 

I previously introduced in our previous work and analyze all the approximations needed to obtain it. 

• , We also provide some extensions of the model to improve its applicability. 



I. INTRODUCTION 

The idea that a spatial soliton and a surface plasmon polariton (SPP) can couple to form the soliton-plasmon, or 
soliplasmon, bound state was originally proposed in Ref. [l| . These authors proposed by means of a sharp physical 
^-j- | intuition that the bound system should obey in certain limit a nonlinearly coupled oscillator model with a peculiar 
nonlinearity, originated by the soliton tail at the metal interface, which was considered as the driving mechanism for 
the coupling and the dynamical features of the soliplasmon system. Despite the amount of undetermined coefficients 
in the model, relevant qualitative predictions of the soliplasmon properties were made. These nonlinear modes where 
later numerically demonstrated to exist in the context of the general vectorial nonlinear Maxwell equations and the 
nonlinear oscillator model (found with asymmetric coupling) was presented with fully determined coefficients @. 
However, that oscillator model was never proved. 

Although the soliplasmon proposal is relatively recent, monochromatic surface waves in the vicinity of the interface 
between a nonlinear dielectric and a metal were extensively studied during the decade of the 80's by many authors 
(see, to cite only a few, pj-[lol|). Amongst those theoretical results, corrections to the profile and wavenumber of i) 
SPP's due to nonlinearity and ii) of the spatial solitons due to the presence of a metallic interface far from its core, 
were presented. At that time, none of the studies, to the best of our knowledge, considered simultaneously a spatial 
soliton and a SPP; the two component wave referred to as soliplasmon. Recently, the model presented in has been 
extended in a first attempt to describe soliplasmons in the challenging and more realistic 2D geometry, consisting 
of several interfaces [Tl|. So far the 2D symmetric like modes are reported, but the antisymmetric ones, which in 
principle require less power in the soliton component 0, remain unfound. On the other hand, the role of nonlinear 
vector effects, generated by strong gradients of the effective nonlinearly induced refractive index, can be of relevance 
for understanding soliplasmon modes near resonance or with high plasmonic component, such as in metals directly 
attached to the nonlinear medium, so standard scalar soliton solutions cannot be appropriate when they are too close 
to the interface and a full- vector soliton solution is required 

In this paper, we present the detailed derivation of the nonlinear oscillator model from first principles. We focus on 
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ID soliplasmons in two geometries of interest, namely the metal/Kerr (MK) 0| and the metal/dielectric/Kerr (MDK) 
[l[ interfaces (by dielectric we implicitly mean a linear dielectric and by a Kerr dielectric with cubic nonlinearity) . 
The paper is organized as follows: in Section II, we motivate and present our variational ansatz and the main features 
of soliton-plasmon coupling in the general context of nonlinear Maxwell's equations; in Section III, we focus on the 
soliton component of the soliplasmon ansatz in the case of total decoupling and we obtain the dynamical equation for 
the soliton variational parameter; in Section IV, we obtain the equation for a general nonlinear plasmon stationary 
mode (decoupled from soliton) and construct the dynamical equation for the plasmon variational parameter. Finally, 
in Section V, we obtain the nonlinear resonator model for the coupled system formed by a nonlinear plasmon and 
a soliton in the weak coupling approximation. In Section VI, we pay attention to the cases of a MK and MDK 
geometries previous mentioned recovering and demonstrating the model presented in Refs.P, 0|- In all sections, we 
stress the role played by the different approximations used in our derivation. The relevance and differences of this 
model with respect other approaches are especially emphasized in the conclusions. 



II. VARIATIONAL ANSATZ FOR NONLINEAR MAXWELL'S EQUATIONS 

We consider the wave equation for the electric component of a monochromatic EM held describing its propagation 
in an optical media characterized by inhomogeneous linear and nonlinear (Kerr) susceptibilities. We take the most 
general form of these equations derived from Maxwell's equations without assuming neither the scalar nor the paraxial 
approximations. Thus, our starting point is: 



V 2 E- V(V-E) = -k 2 T>(E) 



-fc 2 e L E - k 2 P NL (E) 



(1) 



where £l(x) = 1 + x^( x ) is the in-homogenous linear dielectric function and 



P NL (E) = x 16 ' (E • E*) E + X w (E ■ E) E*, 

is the Kerr nonlinear polarization, where %^ 3 ^(x) and x^ 3 H x ) are the in-homogeneous third order susceptibilities 
associated to the Kerr effect. All functions have an implicit, but undisplayed, dependence on the frequency of the 
nonlinear wave uj = cko. 

Our goal is to find a simplified, although physically meaningful, model describing the existence of soliton-plasmon 
resonances, or soliplasmons, obtained by the interaction of a spatial soliton of the Kerr medium with a Surface 
Plasmon Polariton (SPP) on a metal/dielectric interface. We consider a simple configuration consisting in a spatial 
soliton moving (within a Kerr medium) in parallel to a metal/dielectric interface interacting with a SPP propagating 
on it (see Fig[T|). It has been proven that stationary nonlinear states of this system in the form of soliton-plasmon 
resonances indeed exists as a solution of Eq.([T]) for this type of configuration 0|. They show a structure in which the 
soliton and plasmon components are clearly distinguishable, which supports to adopt the following variational ansatz: 



E : 



E„ P {Mz)}Zi 



E, 



(2) 



in which, in principle, we do not explicitly display the dependence on the variational parameters so that a multi- 
parametric dependence can be assumed. We will particularize these parameters later. We also allow the SPP to 
behave nonlinearly by the effect of the Kerr nonlinearity on its own propagation 0, [2l[ thus giving rise to what we 
call a nonlinear plasmon. We introduce the general ansatz in Eq. |1| to get: 



d 2 E np 

dz 2 



L E np - V (V • E np ) 



d 2 V s 

dz 2 



LnEt* 



— — fcgPNL (E np ) — fc 2 PNL (E s ) — fc 2 QK (E np , E s ) , (3) 

where we have taken into account that the soliton is essentially a scalar solution so that V • E s f» 0. We have also 
defined the differential operator Lq as Lq = V 2 + fc 2 £z, where V* = (d x , d y ) is the transverse gradient operator. The 
term Qk (E np , E s ) includes all nonlinear Kerr terms coupling the SPP and soliton components: 

Q K (E np , E.) = Pnl (E np + E.) - Pnl (E np ) - P NL (E.) . (4) 

The SPP and soliton fields have different features. The SPP, as a surface wave and even when it behaves nonlinearly, 
can only exist bounded to the metal/dielectric interface, its nature being purely vectorial. On the other hand, the 
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Figure 1: Parallel illumination of a metal/dielectric interface from a Kerr medium. 



soliton can move freely in the Kerr medium with no restriction and its existence is due to the Kerr nonlinearity and it 
is essentially a scalar wave. It is clear from physical arguments and from results in Ref.fH that the coupling between 
these two entities vanishes as we move the soliton far way from the interface. In such a case, the overlapping of E np 
and E s tends to zero and, thus, Qk — > 0. In the limiting case in which they are infinitely far away, Eq.([T]) gives rise 
to two independent equations for E np and E s : 



d 2 E np 

dz 2 



+ L E np - V (V • E np ) = -fc2p NL (E np ) 



and 



d 2 E s 

dz 2 



L V S = -fc 2 P NL (E.) . 



(5) 



(6) 



Thus we can consider the action of Qk as a perturbation that couples the solutions of these two independent equations. 
For this reason, in our variational approach we will first consider the situation in which the two previous equations are 
satisfied separately. We will write the variational equations for both components separately and then we introduce 
the first order correction originated by the coupling term Qk • 



III. VARIATIONAL EQUATION FOR THE SOLITON 

We start by considering the equation of an uncoupled soliton field satisfying Eq.([6]). We consider the simplest case 
by choosing the fundamental soliton solution as a variational ansatz and by taking a single variational parameter, 
namely, the soliton amplitude. In the planar geometry under consideration, as given by Fig[TJ we consider a x and z 
dependence only, so that the fundamental soliton corresponds to that of the ID Helmholtz equation. This stationary 
solution has the sech form: 

e lf3sZ =uE s {x;C)e l ^ z , C 6 R+ aeR+ , (7) 
where 7 = + > u is a real unitary vector and the propagation constant of the soliton is given by 



E s (a;, z) = uCsech 



UqC (x — a) 



/3 S 2 = kl (e K + \C 2 ), 
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Ek being the linear dielectric constant of the Kerr medium. It is easy to check that the expression ([7]) verifies: 



V? + k 2 (. 



e K +7 E s 



(8) 



or, equivalently, Eq.© for a stationary solution in which Pnl = 7 |E S | Z E S . 

Now, we make the ansatz for the soliton component. As shown in Ref.[2|, for quasi-stationary evolution the dynamics 
of the soliton position a is much slower than that of the amplitude. For this reason, we promote the soliton amplitude 
C to the category of the only variational parameter for the case under consideration. Therefore, we establish the 
ansatz as follows: 



E s (x, z) -> E a (as, z) = uC(z)sech 



/co |C(z)| (a; — a) 



C(z) e C. 



(9) 



Our aim now is to find the dynamical equation for the variational parameter C(z). We immediately recognize 
that the previous ansatz admits the following useful decomposition (taking into account that we can write C(z) — 
\C(z)\exp[i(p s (z)}): 

E s (x,z) = ue i ^E s (x;\C(z)\), 

where E s (x; \C(z)\) is the stationary soliton solution characterized by the real and positive amplitude |C(z)| and, 

., r ; . Since we are interested in quasi-stationary 



1/2 

consequently, by a propagation constant (3 S — ko (ek + %\C (z) | 2 ) 



evolution, we will assume that the dynamics of |C(z)| is much slower than that of the corresponding phase ip s (z), so 
that, d\C\/dz <C dip s /dz. This assumption is supported by numerical simulations 0|. We substitute now the ansatz 
Eq.© into the soliton wave equation ©: 



d 2 E s 

dz 2 



l e s = -kl [( 



which can be written as 



d 2 E s 

dz 2 



x (3) + f (3) 



Vi + H e K + 7|E s 



IE, 



E, 



E s =0. 



(10) 



But, since according to our ansatz, E s (x,z) = u.exp[itp s (z)]E s (x,z) = uE s (x, z), we can see that the differential 
operator in brackets only acts on the stationary solution E s . Since the latter field, on the other hand, satisfies the 
stationary equation (jSJ), this means that the variational field E s also fulfills 



V 2 t +k 2 (e K + 7|E s 



E s 



(11) 



and, therefore, the previous equation becomes: 

d 2 E s (x; C(z)) 
dz 2 



+ P 2 s E s {x-C{z))=Q. 



At this point, we introduce another useful decomposition for the variational field (0): E s {x,z) = C(z)f s (x; \C(z)\), 
where f s — sech [y^^o |C( Z )I ( x ~ a )] ■ This decomposition permits to write the equation above as: 

^ [C(z)f s (x; \C(z)\] + (5 2 s [C(z)f s (x; \C(z)\] = 0. 

In order to obtain a dynamical equation for C(z), we need to project the previous equation with respect to a suitable 
function. A natural choice for the projection function is the stationary soliton solution at the initial propagation 
point, i.e., E s (0) — C(0)f s (x; |C(0)|). Since C(0) is constant, it disappears in the projection process, so that, after 
projection we have: 



dz 2 



C(z) / f s (x,0)f s (x;\C(z)\) 



PI 



C(z) / f s (x,0)f s (x;\C(z)\) 



= 0. 



We introduce now the notation N s = J R f s (0)f s (z) and we immediately recognize that the dependence of A^ s on z comes 
exclusively from its dependence on the modulus of the soliton amplitude |C(z)|. However, due to the quasi-stationary 
approximation d\C\/dz -C dip s /dz, 



d 2 
dz 



— [C(z)N s (\C(z)\)} 



d 2 C(z) 
dz 2 



N s (\C(z)\). 



(12) 
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Therefore, the variational equation for the soliton parameter takes the simple form: 



— C(z) + f3 2 s (\C\)C(z) 



0. 



IV. VARIATIONAL EQUATION FOR THE NONLINEAR PLASMON 

In the case of the nonlinear plasmon component, we follow a similar procedure as for the soliton case. However, 
equations are here more cumbersome to analyze since we have to deal with the vectorial part of the differential operator 
in the wave equation (J5j> - Another difference is that, since the SPP is a surface wave, the linear dielectric function is 
no longer a constant, as for the soliton equation, but rather a function defining the dielectric/metal interface: 



€l{x) = s p (x) 



E m if X < 

£d if x > 0. 



(13) 



A nonlinear plasmon is a stationary solution of Eq.([5|) of the form: 



E np (x, z) = e np (x)e 



i/3 np z _ ( e npt {x) ^ pifs 

^■npz \X*) 



where e npt = (e np:c ,e npa ) stands for the transverse components of the electric field. We will consider that the 
nonlinear plasmon stationary solution is a conservative soliton. This means that for the stationary solution we will 
assume that the system has no losses, neither linear nor nonlinear, so that El will be a real function. For the same 
reason, we will take real nonlinear susceptibilities (x^Kx^ £ R )■ The complex character of a realistic el will be 
taken into account when we set the dynamical equations for the variational solution. 

According to Eq.([5]), the transverse components of the nonlinear plasmon solution verify: 



/3 np e n pt + V e npt 



V t (i/3„ p e npz + V f • e npt ) 



k^EpQ^t 



K 



,0) 



* \ , -(3) / \ * 

^np ' ^npj ®npt T" X V^np ' *^npj G n 



pi 



(14) 



In the linear case, the transverse components of the electric field corresponding to the eigenmodes of an axially- 
invaria nt sy stem can be chosen to be real functions (e npt = e* pt ) whereas the axial ones are pure imaginary (e npz = 
-e* ) [13J. As we will see next, this choice is also consistent in the nonlinear vector case. Assuming these properties 
for the electric field, the transverse component of the nonlinear polarization term associated to the previous equation 
can be written as 



NLi 



x (3) (t> 

,,(3) 



np " e np) e npt 



7l e npt| 



T^npjs 



j npt 



X ( e np " e np) e n pt 

(^np£ ' &npt 



where 7 = x' 3 '' + an d 7 = X ~~ X^ ■ An analogous calculation leads to the following relation for the axial 
component 



-Pnlz = 



7 |e np t| 



7 e np2 



J/3 ap z 



The total displacement vector D = eE + Pnl takes then the form 



D t = 

D z = 



£l + (7|e„ P t| 2 ~ie\ v . 



c npi c — c np c npr c - 



2 2 

ptl -7e nF 



(15) 



On the other hand, due to the mathematical identity: 

V [V 2 E- V(V-E)] =0, 
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it is identically verified from the nonlinear vector wave equation ([T]) that 

VD = 0, 

which imposes a constraint between the axial and transverse components of the stationary solutions: 

i 1 

V t • (£n P e n pt) + i(3 np £n P e np z = => e npz = -r—^V t ■ (e np e np t) ■ (16) 

Pnp ^np 

Despite its form, the previous equation does not provide an explicit expression of the axial component in terms of 
the transverse ones. The reason is the dependence of both e np and e np on e npz as well. In the most general case, the 
solution of the nonlinear vector problem requires to solve the transverse equation (|14[) along with the constraint (|16[) 
in a self-consistent manner. The form of the constraint ([I7i]) also demonstrates the consistency of the assumptions 
with respect the real character of transverse components and the pure imaginary condition for the axial one. Indeed, 
this constraint shows that if e npt E K automatically e npz becomes a pure imaginary function. This is so since we are 
considering el, an d to be real, so that e np and e np also are. 

In many situations, despite the stationary eigenmodes have an hybrid nature, the axial component is commonly 
remarkably smaller than the transverse one. So that, we can reasonably consider in many circumstances that |e npz | <C 
I enpt |- We refer to this condition as the quasi-transverse approximation and, in practice, it will implemented by 
neglecting terms which are second order in the axial component, i.e., 0(e^ pz ) — > 0. 

When the quasi-transverse approximation is considered, the nonlinear vector eigenmode can be described by two 
nonlinear effective functions depending only on transverse components: 

£np ~ £L +7l e npt| 2 , £np ~ £L +7l e npt| 2 - 

In this approximation, transverse and axial components decouple in the equation for e npt Eq. (|14p since, according to 
the constraint (|16l) . it is verified that 

■a ^ n Pw ^t^np /i-7\ 

«Pnpe n p2 ~ — • e npt = • e npt , (17) 

£np £np 

The previous equation permits, after substitution into Eq. (|14p . to eliminate the axial component completely from 
the transverse equation. 

(V 2 + kle nv ) e npt + V t (F npt • e npt ) = p 2 np e npU (18) 

where F npt = e~ p V t e np + 5V t , 8 — (e np — e n p)/£np being the nonlinearly-induced anisotropy function. Note that 
due to the quasi-transverse approximation, e np and F np4 depend on transverse components exclusively. Despite this 
fact, it is important to remark that in this approximation axial components are nonzero. They are simply decoupled 
from the transverse ones. They can be obtained in a simple way from the constraint (|17[) once the problem have been 
solved for the transverse components. Unlike for the general case, the constraint becomes now an explicit expression of 
e n pz as a function of e npt . From the computational point of view, the decoupling of axial and transverse components 
considerably simplifies the calculation of the stationary solution since, then, the simultaneous self-consistent resolution 
of the transverse equation and the constraint can be circumvented. 

Up to now, all the analysis is valid for a general linear dielectric function profile el which does not need to be 
necessarily that of a SPP on a metal/dielectric interface. This means that results can be applied to arbitrary axially- 
invariant structures even in 2D. However, since we are interested in the case of a SPP on a planar structure (FigfT]), we 
will assume that we deal with a ID nonlinear plasmon in a TM configuration. The corresponding electric field has then 
the form e np = (e npx , 0, e npz ) T and thus the transverse vector has no component in the y direction e npt = (e npx , 0) T . 
On the other hand, as in every TM mode, e npx and e npz depend on the x coordinate exclusively. The equation we 
obtained in the quasi-transverse approximation (| 18[) becomes then a single equation for the x component e npx (x) in 
which there is no dependence in the y direction. Remarkably, there are some cases for which there exists an analytical 
solution. That is the case of a planar metal/Kerr structure @- The general form of a solution of the stationary 
transverse problem (fT5)) for a TM mode with only x component has to be analogous to that of the soliton field in the 
previous section: 

E npx (x, z) = e npx (x)e^ z = Af np (x; A)e^ z , A £ M + 

where A is the nonlinear plasmon amplitude. As for the soliton case, we choose it to be the peak value for e npx 
(A = |e np2 ..o|) whereas / np plays the role of the sech function. Inasmuch the transverse equation only depends on e npx 
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and not on e npz we only have a dependence on A and not on the amplitude of the axial component. In the general 
case of a TM mode with coupled axial and transverse components we would have instead: 



-npr 



-npz 



(X) 
(x) 



Af ap (x;A,B) 
Bg np (x;A;B) 



In our case, since we are applying the quasi-transverse approximation, we keep the dependence on A exclusively. It 
will be this coefficient the only one that we will promote to variational parameter. In order to select our variational 
ansatz we proceed analogously as for the soliton held in the previous section. We transform the stationary solution 



E npx into the variational ansatz E npx according to the following rule: 



E npx (x,z) -> E npx (x, z) = A(z)f np {x; \A(z)\) A(z) e 



(19) 



It must be clear now that once the stationary problem has been solved for all values of A (a prerequisite that must be 
fulfilled prior to the analysis of the variational equations) , the function / np in Eq. (|19p is perfectly known for a given 
value of A{z). In some particular cases, such as in Ref.[3] it is even possible to provide an analytical expression for 
the stationary solution and, consequently, also for / np . 

As before, we now introduce the ansatz (|19p into the dynamical nonlinear plasmon equation for the x component 
([5]) to obtain, in the quasi-transverse approximation: 



d 2 E np 
dz 2 



(V? + k 2 e p ) E npx - 



d_ f dE n 
dz \ dx 



p~- 



d_ 

dx 



(V t • E n 



pi) 



-kfa\E npt \ 2 E n 



px 



(20) 



where el = e p is the linear dielectric function profile of the metal/dielectric interface, as in Eq. (fTU|) . We take into 
account now that, according to the variational ansatz ()19j) . we have (we write A{z) as |A(z)| expit/? p (z)) 



E npt (x,z] 



_ e iv P (z) 



\A(z)\f np (x;\A(z)\) 




_ e if> P (z) 



(x;\A(z)\) 




e npt (x; \A(z)\), 



where the function e npt (x; |A(z)|) is the solution of the stationary equation (|18[) with real amplitude |A(z)|. Notice 
that, consequently, E np t(x, z) also satisfies the stationary equation (fT8|) with eigenvalue /3^ p (|>l(z)|). The function 
e np t verifies the constraint (|17[) in the quasi-transverse approximation. Thus, due to our variational ansatz, we can 
also find the corresponding constraint for for the variational field E npt : 



V,-E. 



(21) 



Introducing the constraint above in Eq. ([20|) . we obtain 



9 -Enpx 

dz 2 



d_ 

dz 



dE np2 
dx 



dE npz 
dx 



+ [V 2 t + k%£ np ] E npx 



(F np t ■ E npt ) 



(22) 



Taking into account that E npt (x, z) verifies the stationary equation (TT5)) with eigenvalue /3^ p (|>l(z)|), the equation for 
the variational field experiments a notable simplification 



d 2 E 
dz 2 



d_ 

dz 



dE n 



dx 



. R dE npz 
iPnp dx 



PlpE n 



0. 



(23) 



Now we can proceed to project this equation onto a suitable function in order to find the dynamical equation for 
the variational parameter A(z). Since we are dealing with a vector equation for the electromagnetic field, we should 
perform the projection using a proper scalar product for this case. Orthogonality for vector eigenmodes is defined 
through the vector relation J (H* x E) ■ z which equals J H y E x in our case. Thus, we should project Eq. ([2"3")) onto a 
suitable selected value of H y (x,z). As for the soliton case, a natural choice is to take the stationary solution at the 
initial propagation point H y (x,0). After performing the projection, we find 



dz 2 



H y {0)E npx (z) 



~ I Tz ~ l/?np 



Hy(0) 



dE n 



dx 



2 



H y {0)E npx (z) 



= 0. 



First of all, we pay attention to the integral involving the axial component: 



d 

H y (x, 0)— E npz (x, z) 



H y (x,0) E npz (x,z) 



d_ 



i— \ e np E npz (x,0)E npz (x,z) 
c Jr 

0(E 2 npz ) -> 0, 



(24) 



which vanishes in the quasi-transverse approximation. We have used here the Maxwell's equation [V x H] = 
dHy/dx — —iu)c~ 2 eE z to write the previous integral in terms of the axial component of the electric field. 

Now we consider the realistic situation in which the system admit losses, a circumstance which is immediate in 
the case of metals. We return to the dynamical equation for the variational field (|2"2"|) and consider now that e p is a 
complex function, so that, we make the substitution e p — > e p + iAie p , where we keep in our notation e p as the real 
part of the dielectric function whereas A;e p is a function indicating the distribution of linear losses in the system. 
This substitution generates an extra term in Eqs. (|22[) and (|23|) of the form ik$ Ai£ p E npx , in such a way that, after 
the projection, we obtain: 



dz 2 



H y (0)E npx {z) 



■B 2 



H y (0)E npx {z) 



I 

J J: 



H y {Q)&ie p E npx {z) 



0. 



We recall that for a TM stationary mode it is true that 

k a c 



Hy(x,0) 



E npx (x,0) =Kf np (x; L4(0)|), 



where if is a constant independent of x and z, so that it disappears from the equation. On the other hand, 

E npx (x, z) = A(z)f np (x; \A(z)\). 

Therefore, 



where we have defined the "norm" N np as 



— [N np A(z)] + N np [B 2 ip + i\ Dp ] A(z) = 0, 



iVn 



and the loss parameter A np as 



Ann = 



/ np (a:;|A(0)|)f n p(a:;|A(«)|) 



f np (x; |A(0)|)A Z £ p (a;)/ np (a;; \A(z)\). 



(25) 



(26) 



By construction, the "norm" 7V np and the loss parameter A np depend on the modulus of the variational parameter 
|^4(z)|. Likewise, the propagation constant B np of the nonlinear plasmon is also dependent on this quantity. All these 
dependence on |A(z)| can be fully established once the stationary nonlinear problem has been thoroughly solved. 
A further simplification of the variational equation is permitted in the quasi-stationary approximation d|^4|/rfz <C 
difp/dz, since we can proceed as we did for the soliton case to write 



— [N np (\A(z)\)A(z)]^N np (\A(z)\ 



dz 



so that we finally obtain 



r ^2 



A^np 



—A(z) + B 2 (\A\)A(z) 



d 2 A(z) 
dz 2 



= 0, 



where we have defined a complex effective nonlinear propagation constant as f3 2 p = B 2 p + ?A np . 



9 



V. VARIATIONAL EQUATIONS FOR THE SOLIPLASMON BOUND STATE 

In this section we will establish the variational equations describing the propagation of a soliplasmon bound state. 
Our variational ansatz for a soliplasmon will be given by a superposition of a nonlinear plasmon and a soliton. 
However, instead of assuming a general dependence in multiple variational parameters, as in Eq.@, we will reduce 
the number of variational parameters just to two: one associated to the nonlinear plasmon, A(z), and a second one, 
associated to the soliton, C(z). They correspond to the amplitudes of the nonlinear plasmon and soliton solutions as 
defined in Eqs.© and ([19]), Therefore, our variational ansatz for the soliplasmon solution will be: 

E(x, z) = E np (x; A (z)) + uE s (x; C (z)) . (27) 

In the same way, we will treat mathematically the plasmon and soliton components as we did in the previous two 
sections. This means that we will assume the same approximations we used to demonstrate the variational equations 
for the uncoupled system. Summarizing, we will work under the following approximations for the variational fields: 

• Scalar approximation for the soliton field: V • E s w 0. 

• Quasi-transverse approximation for the plasmon field. Terms of order E 2 pz and higher will be neglected: 
O(E 2 npz )^0. 

• Quasi- stationary approximation for both: d\C\/dz -C d(p s /dz and d\A\/dz <C d(p p /dz. 

A. Equations for the coupled plasmon and soliton variational fields 

Taking all these approximations into account we proceed to find the corresponding variational equations for the 
plasmon and soliton field components of our ansatz (1271) . As mentioned in Section UH when the soliton is located 
infinitely far away from the nonlinear plasmon, the two localized solutions in Eq. (|27[) present a vanishing overlapping, 
so they can be treated independently in such a way they verify uncoupled independent equations. This is the analysis 
we have performed in the two previous sections. When this overlap cannot be neglected, an explicit coupling appears 
and then Eq.© holds instead. On the other hand, the soliplasmon solution found in Ref.jH is a TM mode of the 
electromagnetic field. Besides, the axial component of its electric field is only relevant for the plasmon field close to 
the interface and not for the soliton. For these reasons, we consider that E = (E x ,0,E z ) and u « (1,0,0). In this 
way, the axial component will be approximately given by the plasmonic component exclusively (E z w E npz ) so that it 
will be possible to evaluate it from E npx according to the procedure presented in the previous section. Consequently, 
our starting point will be the variational equation for the x component of the electric field. So that, we write Eq.Q 
for the x component (recall that Lq = V 2 + Ii^el) 

d 2 £n P , d fdE npz \ dE Dpz d d 2 E 

~d^~ " Tz [~^x~ ) + lP ^~dx~ + LoEnpx + d~x (Fnpt ' Enpt) + ~dz^ + L ° Es ~ 

-k 2 7 |E npt | 2 E npx - k 2 al \E S \ 2 E s - k 2 Q K (E npx , E.), (28) 



where we have used the constraint ([21]) for V • E np — valid in the quasi-transverse approximation — and defined 

Qk = Qk,i- 

There are two coupling mechanisms in Eq. (|28p . One is purely nonlinear and it is generated by the Kerr coupling 
term Qk- The other one is due to the presence of a variation of the linear dielectric function in the regions where the 
field is localized, i.e, in the nonlinear plasmon and in the soliton regions, in our case. This is a well known mechanism 
in solid state physics and it is the origin of the coupling between neighboring wave functions in the so-called tight 
binding approximation [14j |. Let us see how it works in the present case. We note that the linear operator Lq in 
Eq. (|28p does not coincides exactly with that corresponding to the uncoupled solutions in the variational equations for 
the soliton (Eq. (fT0|) ) and nonlinear plasmon (Eq. ([22|) ). The reason is that the total linear dielectric function differs 
from the uncoupled ones in the regions of the ID space where the functions are not localized. To be more specific, let 
us write this function for a MDK structure as in FigJTJ 



e L ( x ) = }W (29) 
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where d is the thickness of the dielectric layer. In this case, e p represents the profile of the linear dielectric function 
for the plasmon component, i.e., it defines the dielectric constant profile of the MD structure, as defined in Eq. (|13p . 
We could consider, with any lack of generality, that there is also a modulation of the linear dielectric function in the 
Kerr medium, so that, we would have a function e s (x) instead of £k in the previous expression. However, in order to 
preserve the analysis of the soliton component exactly as we did in Section [HII we will keep ek as the linear dielectric 
function for the soliton field. Generalization to arbitrary e s (x) will be straightforward once final results are obtained. 
The definition of the total dielectric function in Eq. (f2T))) suggests the following two decompositions for el- 



e p + Ae p 
e K + Ae s 



(30) 



where the local variations Ae p and As s would be given by 



As p (x) 



x < d 

£k — £p(x) x > d 



x < d 

£k — £d x > d 



(31) 



and 



As s (x) 



£ P (x) — ck x < d 
x > d 



e m - s K x < 
£d — £k < x < d , 



(32) 







x 



> d 



where we have taken into account the dielectric function profile for the MD interface as given by Eq. (fT3"|) . The 
decomposition (|30[) of the linear dielectric function permits to write the operator Lq in two different ways: 



Lo = (V? 



k^Ep) 



L = (\7 2 t +kls K 



Ae p — Lq p 

- Ae s = L 0s 



Asp 
-As s 



and, therefore, we can rewrite Eq. (|28l) as 



d 2 E n 



\).r 



dz 2 



d_ f dE npz 
dz \ dx 



+ ■ 



d 2 E s 
dz 2 



dE n 



pz 



Ox 



L 0s + k 2 j\E s 



E, 



cigcvaluc equation for soliton 



J 0p 



2l 9 
|E np t| E npx + — (F npt • E npt ) 



cigcvaluc equation for NL plasmon 

-kl {Ae p E npx + Ae s E s ) - k 2 Q K (E npx , E s ). 



We immediately recognize in the previous expression the appearance of the nonlinear operators for the soliton and 
plasmon stationary solutions. We demonstrated that soliton and plasmon variational fields were also eigefunctions of 
these operators, so that a simplification of the above equation can be obtained: 



d 2 E npx ( d 



8z 2 



dz i/3np 



dE n 



d 2 E s 



-kl {Ae p E npx + Ae s E s ) - klQ K (E npx , E s ). 



(33) 



B. Dynamical equations for the variational parameters in the weak coupling approximation 

In order to obtain the equations for the variational parameter A(z) and C(z) we need to project out Eq. (|33[) into the 
adequate projection functions. In the two previous sections we projected the two uncoupled equations for the soliton 
and plasmon variational fields using suitable soliton and plasmon projection functions for each case. Here, we will 
make an identical choice for the projecting functions and we will make two different projections of Eq.(|33j). The first 
one, corresponding to the soliton projection, will be performed with respect to the soliton stationary field at z = 0, 
i.e., E s {x, 0). For the second one, corresponding to the plasmonic projection, we will use the plasmon stationary field 
at z = 0, i.e., H y (x,0). 

When performing the aforementioned projections in Eq. (1331) , we will encountered a new type of overlapping integrals 
not present in our previous analysis. They correspond to integrals involving products of soliton and plasmon fields. 
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Since soliton and plasmon functions are localized in different regions of the space, these integrals are expected to be 
small when the soliton field is localized sufficiently far way from the MD interface. More specifically, if we analyze an 
overlapping integral of the form (we assume F p to be a plasmonic function tightly localized around the interface) 

In(a)= [ F p (x)f?{x-a), 
Jr 

and we consider the overlapping to be small (a 3> 1, implying |7„| <C 1), we can estimate the order of this integral 

1/2 

by approximating the sech function by its exponential tail, so that, f s w exp« s (x — a)), where k s = (7/2) fco \C\. 
Therefore, 

I n (a) » e~ nK ° a [ " F p (x)e nK » x ~ 0{e~ K " a ) n , (34) 
Jo 

where d p is a small parameter of the order of the penetration length of the localized function F p into the dielectric 
medium. The previous argument suggests to take e s = e - K " a as a small parameter. Analogously, terms depending on 
the plasmonic tail of the form e~ Kpd where the distance d is substantially larger than the plasmon penetration length 
(d 3> «p 1 ) will be also small, which defines e p = e~ Kpd as a small parameter as well. Hence, we define the following 
additional approximation associated to the coupling of soliton and plasmonic components: 

• Weak coupling approximation. 

— Terms of order e~ 2Ksa and higher will be neglected: 0{e~ 2KBa ) = 0{t 2 ) — > 0. 

— Terms of order e~ 2Kpd or higher will be likewise neglected: 0(e~ 2Kpd ) = 0(e 2 ) — > 0. 

The previous approximation complete the set stated previously. Now, we use it to perform the projections by keeping 
only the leading terms. 

We can further simplify the soliplasmon equation for the variational fields (|33[) by simultaneously invoking the 
quasi-transverse and weak coupling approximations. Let us pay attention to the second term depending on the axial 
component E npz . The plasmonic projection of this term provides the overlapping integral already seen in Eq. (|24D . 
which, it was proven to be 0(E 2 pz ) and, thus, negligible. The soliton projection provides, in turn, the integral 

I E npz (x)f s (x - a) « e~ K - a [ ' E npz (x)e nK ° x ~ 0(e~^ a )0(E npz ) -> 0, 
Jr Jo 

which we also neglect as a product of two infinitesimals. Therefore, this term does not contribute to any of the two 
projections and it can be also neglected. 

However, this is not the only simplification that we can perform using the weak coupling approximation. They also 
apply to the nonlinear coupling term Q\^{E npyi , E s ). From its definition Q we have 

Qk(£„ p *, E.) = 7 [2\E s \ 2 E npx + E 2 E* npx + 2\E npx \ 2 E s + E 2 npx E* s ] . (35) 

We immediately see from Eq. (|34j) that both the plasmonic and soliton projections of the two quadratic terms in E s 
of Qk are, at least, 0(e~ 2Ks °") [29(. Therefore, the first two terms in the previous equation can also be neglected. 

Taking into account all the previous approximations, our final equation for the variational fields takes a relatively 
simple form: 

r) 2 K r) 2 K 

Q z 2 + Pnp^npx T ^ 2 "t" P s &s ~ 

-k 2 {Ae p E npx + Ae s E s ) - k 2 Ql [2\E npx \ 2 E s + E* px E* s ] . (36) 



Our final step is to project Eq. (|36l) into the soliton and plasmon projection functions in order to convert this equation 
in two different equations for A(z) and C{z). 

1. Plasmon projection 

We start now by performing the plasmonic projection first. As in Section HVl we project with respect to H y (x, 0) = 
Kf np (x; |A(0)|), where the constant K can be ignored since it disappears after the projection. We obtain 



A^np 



3ps 



= -A PP A - A ps C - A K [2\A\ 2 C + A 2 C*] , (37) 
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where N np is defined as in Eq. (|26l) and we have defined three new coupling terms 

S ps = f /„ P (0)/ S (z) 
A ps = kl £f np (0)AeJ s (z) 

A K = k 2 [ 1 f np (Q)fl p (z)f s {z). (38) 

JR 

However, the plasmonic self-interaction term A pp term can be neglected according to the weak coupling approximation 
since 

A pp = kl f f np (0)Ae p f np (z) ~ k 2 (e k - e d ) f e~ 2K - x ~ 0{e~ 2K ? d ) -> 0, (39) 

JR J d 

where we have taken into account the form of the local variation of the dielectric function for the plasmon as in 
Eq. (|3ip and we assume that the width of the dielectric layer d is larger than the plasmon penetration length in the 
dielectric: d ;§> k^ 1 . Besides, we have approximated / np by its linear counterpart / np « f p ~ e~ KpX because in the 
regions where de plasmon field is very weak nonlinear effects are negligible. 

In order to obtain the plasmonic projection (|37[) we have also used the quasi-stationary approximation for the 
plasmon field d\A\/dz <C dip p /dz that we already used in Section llVl to write d 2 (N np A)/dz 2 « N np d 2 A/dz 2 . An 
analogous argument, in this case using the quasi-stationary approximation for the soliton field d\C\/dz dip s /dz, 
has been also used to prove that d 2 (5 ps C)/dz 2 w S ps d 2 C/dz 2 . 



2. Soliton projection 

In this case, we project Eq. tptr?)) with respect to the stationary soliton field at z = 0. This field is given by the same 
projection function, E s (0) = C(0)/ s (0)), that we used in Section Hill Since C(0) is a constant not depending on x nor 
z, it disappears after the projection is carried out, so only the spatial function f s (x; |C(0)|) appears in the projection 
integrals. In this way, the resulting equation obtained after soliton projection is: 

= -A SS C - A sp A - A' K [2\A\ 2 C + A 2 C*] , (40) 

where slightly different coupling terms from those appearing in the plasmonic projection — Eqs. (|38[) — are obtained: 

S sp = [ f s (0)f np (z) 

JR 

A sp = k 2 Q f / s (0)Ae p / np (z) 

JR 

A' K = k 2 f 7 / s (0)/ n 2 p (z)/ s (z). (41) 

JR 

The linear self-coupling term A ss can be approximated in the weak coupling approximation as follows 

A ss ee k 2 [ f s (Q)Ae s f s (z) ~ k 2 [" [ Ep (x) - e K ] e K(o)+^](x- a ) _ 0(e -[«.(o)+«.]aj ^ 0; 

JR J -oo 

and, therefore, neglected as its plasmonic counterpart. In the previous integral we have approximated the sech function 
by its exponential tail. This is justified because the integral only covers the metal and dielectric part, so that if the 
soliton field is located not too close to the dielectric region, its value in the integral domain will be given by its 
exponential tail. Using identical argument, it is proven that 

- o( e -K(o)+«.]«) _j. 

Obviously, it is implicitly assumed that the varying value of k s (z) always satisfies the weak coupling condition. So, if 
the condition is satisfied by the field at z = — determined by the value of k s (0) — , then k s (z) ~ n s (0) for all values 
of z so that: 

0( e - K ^ z)a ) - 0(e-^°) Q ) V«, 



d? A + ^ A 



dz 



;C + p 2 s C 
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and, therefore, 

O( e -M0)+K.]a) _ O ^-2«.(0)oj _ 0(e -2 "' ) ^ Q 

Consequently, the weak coupling condition for the soliton field in Section IV Bl should be understood in the above 
sense. 

As for the plasmonic projection, we have also used the quasi-stationary approximation to write d 2 (S sp A) / dz 2 w 
S sp d 2 A/dz 2 and d 2 {N s C)/dz 2 « N s d 2 C/dz 2 . 



3. Equations for the variational parameters A(z) and C(z) 



The plasmonic and soliton projections above provide us already with dynamical equations for the variational 
parameters, but they are not yet in their "canonical form." In order to see this feature, we write them again together 
using matrix notation 



S, n N„ 



A 



App ^\ps 




'A' 




^sp ^ss 




C 





A K AC* 2A K \A\ 2 
A' K AC* 2A' K \A\ 2 





'A' 




C 



Although we know that some of the coefficients in the previous equation vanish in the weak coupling approximation, 
as we have just proven, we will keep them in order to analyze some interesting particular case in the next section. So, 
by pre-multiplying the equation above by 



A n p &ps 
S,n N. 



we can set the variational equations in their standard form: 

Bl p )A =- qp S C-q K (2\A\ 2 C + A 2 C*) 
+ B 2 )jC =- q sp A-q^(2\A\ 2 C + A 2 C*). 
We can recognize that there exist three type of terms in the previous variational equations: 



dz 2 

#_ 

dz 2 



(42) 



Terms modifying the phase velocity of the plasmon and soliton components by renormalizing their propagation 
constant through terms lineal in A and C, respectively: 



-^np 


= P 2 - 

— rnp 


f 3ps^ps 


- N S A PP 


V ^npA s 


0~ps$sp 


B 2 S 




A - 

u sp L -*sp 


- N np A ss 


K N np N s 


0~ps$sp 



• Terms coupling plasmon and soliton components, which are linear in A and C, respectively: 

'A, A, 



*ps 



8ps ^-ss 



N np N s - 8 ps S sp 



Qsp 



$ sp^pp 



N np N s - SpsSsp 

Terms coupling plasmon and soliton components, which are nonlinear in A and C, respectively: 



(43) 



<7K 



N S A K - 


-S PS A' K \ 


N np N s 


3ps$sp J 


N np A' K 


- S sp A K \ 


N np N s 


— Sp S 8 S p J 
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Thus the evolution equations for the variational parameters of the soliplasmon problem (|42[) can be written as a 
system of coupled nonlinear resonators with special characteristics. The special features we are referring to have to 
do with the nonlinear dependence of all the coefficients in the previous equations on the modulus of the variational 
parameters \A\ and |C|. All coefficients depend on quantities which are given in terms of integrals of the stationary 
functions f np (x;\A(z)\) and f s (x;\C(z)\). Note that these nonlinear terms are different from those owned by the 
nonlinear plasmon and soliton when they are decoupled. In the absence of coupling, the nonlinearities come from the 
dependence of /3 np and f3 s on \A\ and \C\, respectively (as we have seen in Sections Mil and IIV1 ) Note as well that 
there are nonlinear dependences which are not explicitly given by the obvious crossed Kerr terms in Eqs. (|42l) . For 
example, the linear coupling coefficients q ps and q sp present nonlinearities which are not directly related to the Kerr 
coupling but to the overlapping of the localized nonlinear plasmon and solution functions / np and f s . 

In order to keep our discussion in the most general form, we have retained terms that do not appear in the weak 
coupling approximation, at least to leading order. If we keep just leading terms in this approximation, the previous 
coefficients considerably simplify. Let us briefly recall the vanishing terms in this approximation: 

• S ps S sp , Sp S A sp , S S pAp S , 8 ps Ak, SspAx ~ 0(e~ 2Ksa ) 0. This is so because all coefficients involved in these 
quadratic products correspond to overlapping functions of /„, and, therefore are 0(e~ Ksa ). 

• Besides, as proven before, A ss ~ 0(e _2Ksa ) — > and A pp ~ 0(e~ 2Kl ^ d ) -} 0. 

• And, finally, also proven before, ~ 0(e~ 2Ksa ) — > 0. 
With all these approximations in mind we find that: 

np ~ r^np 

Bl « Pi 



and 



as well as 

5K 



^ps 
iVnp" 

~n7' 



Ak 

N np 

Qk « 0. 

Summarizing, the general form of the dynamic equations for the soliplasmon variational parameters in the leading 
order of the weak coupling approximation is given by 



VI. CASES OF INTEREST 



The equations found for the evolution of the variational parameters A(z) and C(z) are, in principle, not restricted 
to the specific linear profiles of a MDK structure. Despite we have used the specific form of el for a MDK structure 
in some previous steps to justify some approximations, this procedure has been adopted more for clarifying purposes 
than for necessity. In fact, it is not difficult to realize that the form of the equations and coefficients is preserved if 
we assume two general, localized, linear dielectric function profiles for the plasmon — s p (x) — and soliton — e s (x) — 
fields. 

However, since MDK and MK structures are the main physical motivation of the current study, in this section, we 
will particularize the general approach to these two cases of interest. We will work, in principle, at leading order of 
the weak coupling approximation, so that Eg. (1441) will be our starting point in this section. 
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A. Coupling of a linear plasmon and a soliton in a MDK structure 



We will assume here that we work with a MDK structure, as presented in FigJTl characterized mathematically by 
the function El{x) as presented in Section HVl Eg . ([!?£))) . The plasmon dielectric function e p (x) will be that of a MD 
structure, as given by Eg. (1131) . On the other hand, we will assume that Kerr nonlinearities do not affect the plasmon 
component, something that can be realistically realized by suitable playing with the width of the dielectric layer d 
together with the amplitude of the SPP field. Therefore, considering the plasmon component to be linear implies that 
(since \A\ < 1) 



B 2 kB 2 

as well as suppressing the nonlinear coupling term 0(A 2 ) in Eg. ([4"4")) because 



Hence, we get 



Ak~7/ / P (0)/ p 2 (z)/ s (z)~ 7 



U^ 2 



ipX j 



0{e 



-> 0. 



N, 



-C 



dz 2 



Ps)C =- 



(45) 



It is possible to obtain explicit expressions for A ps and A sp in the MDK case. This is possible because the functions f p 
and f s are known explicitly. According to the definition given in Section [IV! the plasmonic function f p is normalized 
by the peak value of the x component of the linear SPP electric field E px . Since we are dealing with the linear 
solution, we have an analytical expression for it [HI, [l(| : 



E px (x) 



x<0 



PpEg K m x 

PpEo e - Kd x X>Q 
. KnSd 



The peak value of the SPP field is achieved at x — 0, so that 



fp( x ) 



E px (0) 



e K m X x < 
e -K d X x > Q_ 



1/2 

Analogously, the f s function is given by f s (x) — sech [k s (x — a)}, where n s — (j/2) fco|C|. We can approximate the 
sech function by its exponential tail when the overlapping with the f p function is small, hypothesis which is justified 
in the weak coupling approximation. So that, in the overlapping region 

/ S (i)«2e- K ' l " ) , x<d, a>d>0. 



Now, according to their definitions in Eqs (|4"Tj) and ([55)) and to the form of the local variations of the linear dielectric 
functions Ae p and Ae s in Eqs. (f3"Tj) and (|3^|) . we obtain for the coupling coefficient 



/0 pd 
fpfs + {e d ~£K) / fpfs 
-oo JO 



£i<g|Kd — «s 



2e" 



so, 1,2 



+ {£d - £K)d 



(46) 



where in the last step we show an expression valid for sufficiently thin dielectric slabs. 



For A sp we have 



k 2 (e K - £d) 



2k 2 . 



Kd 1^3 



I sip 
^—dKd + {d—a}K s 



(47) 
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where in order to give the approximated expression above we had to assume that > k s and a > rf, so we could 
approximate / s by its right-hand-side exponential tail. 



In the same way, the plasmonic linear "norm" N p = J R f£ can be evaluated to give 



A<„ 



1 



whereas the solitonic one N s = j R f£ takes the simple expression 

2 

Ks ' 



(48) 



(49) 



All the evolution equations we have obtained up to now are second order in the derivative with respect to the 
propagation variable z. Physically speaking, they are non-par axial. However, the physical configuration under 
consideration, in which the soliton propagates in parallel to the metal/dielectric interface (see FigO]) following the z 
axis, exhibits a clear paraxial character. This fact indicates that a slowly varying approximation for the plasmon and 
soliton components with respect the propagation parameter z is expected to be adequate for the analysis of this case 
as well. Certainly, it will properly describe propagation in the regions where most of the energy is localized, namely, 
those where the plasmon and soliton components evolve. Since we are working in the weak coupling approximation 
the plasmon and soliton regions are, by construction, clearly distinguishable in the form of weakly coupled plasmon 
and soliton modes propagating along the z axis. However, despite these components are essentially paraxial, the 
energy exchange between them induced by the coupling is not necessarily so. It can include, in prin ciple, non-paraxial 
components associated to a very fast, in the sense of rapidly varying in z, energy exchange [30]. Thus, we expect 
that, with the exception of this type of rapid oscillations, the slowly varying approximation will provide also a good 
approximation to the solution. Thus, we introduce the slowly varying plasmonic A(z) and soliton C{z) envelops in 
the usual way, taking «k = £f/ 2 as ^ ne reference index, 

A(z) = A(z)e tkonKZ 
C(z) = C*(z)e lfc °" KZ , 

so that, we can write (after neglecting d 2 A/dz 2 <C ik ni~ dA/dz , idem for C), 

—i— = UpA + qC 
dC 



dz 



Li s C + qA, 



in which we have defined the paraxial propagation constants /j, p and /i s as 

_ (P 2 p -k 2 Q e K ) 



2k n K 
(/3 2 s - k 2 e K ) _ fe 7 \n 



2feo^K 4riK 

In the last equation we have used the explicit expression for the Helmholtz soliton propagation constant /3 S 

&o (^k + 2 | C | 2 ) , as introduced in Section Hill 

Analogously, we have defined the paraxial coupling coefficients q and q as 



fc 



2koTiKN p 2ukN p 



fc 



2k n K N s 2n K N s 



fpAe s f s 



fs^^pfp- 



Explicit expressions for the paraxial couplings can be given in the same way as for their non-paraxial counterparts in 
the case of a MDK structure. From expressions (j46|) and (|47| for A ps and A sp , we have 



i%kN p 

fc 

1 

n K N s 



(e m — ess) 

(£K - Ed) 



1 



l 

Kd Kg 



£d — £K 

e m — £k 



-d 
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An important consequence of the previous analysis is the fact that the plasmon-soliton coupling is asymmetric since, 
in general, q ^ q. The previous expressions allow us to obtain relevant information about the characteristics of the 
coupling ratio q/q. At this point, we are only interested in estimating the order of magnitude of this ratio, so that if 
we consider only the leading terms in d, the explicit expressions for N p and N Sl and the relation between the inverse 
penetration lengths in the metal and dielectric, K m and Kd, and the MD dielectric constants [16], 



-koe m \j ■ w k (-£ m ) 1/2 

£m + £d 



Krf = k s d \ — « k e d (-£ m ) 1/2 , 



£m + £d 



we can simplify the ratio into 



q k s K m (e d — e K ) k s ( e d — £k 



q 4n d K d (s m ) K d \ Ae d 



(50) 



where we have simultaneously assumed that K m 3> n d 3> k s and \s m \ 3> £ dl £k- The latter approximation is 
pretty realistic for common MD interfaces. The former one simply indicates that the plasmon penetration length in 
the metal is significantly smaller than in the dielectric (which is consistent with the previous approximation since 
Km/ftd — £m/£d 3> 1) whereas the plasmon penetration length in the dielectric is, in turn, smaller than the typical 
spatial soliton width (which is also a reasonable assumption for paraxial solitons). 

The previous result shows that this ratio is generally small by two reasons: first, and most important, in this regime 
K s /n d <C 1, and, second, for standard dielectric and Kerr materials the dielectric constants ratio (term in parentheses 
in Eq. H50p ) can be also small. Let us be more precise and introduce explicit expressions for n d and n s . For the soliton 

1/2 

inverse penetration length k s = (j/2) 1 fco|Co| we choose the peak amplitude to be the initial one since, in order 
to preserve the weak coupling approximation, it must be true that K a (z) ~ K s (0) for all values of z. It is useful to 
introduce the dimensionless nonlinear coefficient 7 = -7 1 Co | 2 , normalized to the initial soliton peak amplitude. In this 
way, we get 

J/3 

q / £d - £K \ / -£m7 



q V 4 £ 2 

If we introduce some typical numbers for the dielectric constants, say e d = 1.5 2 , £k = 2 2 and e rn = —80, we can 
make the following estimation 

I ~ 10- 1 (407) 1/2 ~ 0.57 1/2 . 

q 

This ratio is small because in ordinary cases 7CI. Now, for a standard nonlinear Kerr medium with a nonlinear 
index n-i ~ 10 _19 (m 2 /W) (one order of magnitude higher than that of silica), corresponding to a value of 7 = 
ceo£K^2~10 _22 (m 2 /V 2 ), and a peak electric field Co = Eq ~ 10 7 (V/m), which would be a typical value for a 10 fim 
wide (along x direction), 100 /^m (or larger) high beam (along y direction) — as an approximation to a ID soliton — , 
with peak power of Po > 10 kW, we would get 7 ~ 10 -5 . The order of magnitude of the coupling ratio would be then 
q/q ~ 10 -3 which shows how small this ratio can be in typical situations. 

In summary, the paraxial equations for a soliplasmon state in a MDK system in which the plasmon component 
behaves linearly are, under all the approximations carefully explained in the present section and in matrix form, 



dz \C J \ q J \ C 

in which, in typical conditions, the soliton-to-plasmon coupling q is much stronger that the plasmon-to-soliton one q: 
q <C q. This equation is the one presented in Ref.@ in which c p = A and c s = C are the paraxial plasmon and soliton 
envelopes, respectively. 



B. Coupling of a nonlinear plasmon and a soliton in a MK structure 



We analyze now another interesting case, namely, that of a metal surface directly attached to the Kerr medium and 
subject to parallel illumination by means of a spatial soliton of the Kerr medium (see Fig(2]). In this situation, the 
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Figure 2: Parallel illumination of a metal/Kerr interface from the Kerr medium. 



dielectric slab of the MDK structure analyzed in the previous section is no longer present. This has two important 
consequences in the mathematical description of the problem: firstly, the linear index contrast between the dielectric 
and Kerr medium of the MDK structure disappears. e e i — Ek] and, secondly, because the Kerr medium is now attached 
to the metal, Kerr nonlinearities directly affect the SPP, so that it is required to consider a nonlinear plasmon instead 
of its linear counterpart. Let us see how these two features affect the soliplasmon propagation equation. 

As before, we work at leading order of the weak coupling approximation, so Eq. fB)) is the correct variational 
propagation equation to use in this case. The first consequence of working with a MK structure instead of with a 
MDK one is that the plasmon-to-soliton coupling completely disappears. Indeed, according to Eq. ([3T]) Ae p — and, 
therefore, A sp — according to the definition (j4Tj). On the other hand, the plasmon propagation constant /3 np is now 
that of a stationary nonlinear plasmon and it depends on the SPP amplitude The analytic form of /3 np on \A\ is 
know only in some cases [3j. In the general case, /3 2 p is nothing but the eigenvalue of the nonlinear equation for the 
electric components of the nonlinear SPP field (fl~8|) . As explained in Section llVl the nonlinear terms in this equation 
depend nonlinear ly on the plasmon amplitude as \A\ 2 . So that, /3 2 p — (3 2 p (\A\ 2 ) and we can always perform a Taylor 
expansion in \A\ 2 : 



We are interested here in the first order nonlinear corrections to the linear case, so that we keep the first correction 
only and neglect 0(|A| 4 ) terms. Consequently, to leading order in the weak coupling approximation and to first order 



+ 0(|A| 4 ). 



in \A\ 2 , we have 




(51) 



We can write the previous equation in a form that resemble that of a linear plasmon in Eq. (|45p . 





^+(3 2 (\C\)]c = 0, 



(52) 



where now 




AC* 



(53) 
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Note that the effective propagation constant is now a complex number /3 np S C because of the nonlinear AC* coupling. 
Thus the nonlinear soliton-to-plasmon coupling simultaneously affects both the phase velocity of the nonlinear plasmon 
and the plasmon amplitude, due to the presence of a nonzero gain-loss coefficient proportional to the imaginary part 
of the effective propagation constant: 3(/3^ p ) ~ sin(ip p — ip s ) ^ 0, in general. There is also a nonlinear plasmonic 
modification of the coupling coefficient A ps in Eq. (|53p that can slightly modify the nature of the coupling with respect 
the pure MDK case. These two effects are additional to those associated to the ^-independent coefficient A ps (|C|) 
appearing in the case of the linear-plasmon/soliton coupling in a MDK structure [2J. 

We see that, at this order of the weak coupling approximation, the soliton equation decouples from the plasmon one. 
Physically speaking, the soliton acts as a non-depleting reservoir pumping the non-linear plasmon without experience 
any energy exchange with the SPP. Obviously, this is only true within the order of our approximation because we 
have truncated higher-order terms. Indeed, there exists this type of energy exchange from the plasmon to the soliton 
but this has to fulfill two conditions: (i), it has to be very small, i.e., in mathematical terms is has to be at least 
0(e~ Ksa ) 2 since this is the order of terms neglected in our weak coupling approximation; and (ii), it necessarily has 
to arise from terms neglected in the process of deriving the variational equations at leading order, since the leading 
order plasmon-to-soliton coupling for a MK structure, linked to variations in the linear dielectric functions, is strictly 
zero inasmuch Ae p — 0. So, we would need to go to next-to-leading order in the weak coupling approximation if we 
wanted to account for these effects. 

This generalization to next-to-leading order is certainly more complicated because we need to include terms that 
we have neglected in previous analysis. If we keep all the approximations to the same order, as before, and we only 
go to the next order in the weak coupling approximation, we need to consider the four terms in the Kerr coupling Qk 
defined in Eq. (P5j) . We have to retained now the first two terms. We neglected them before because they provided 
terms 0(e~ Ksa ) 2 that we want to keep now. The plasmon and soliton projections provide two more terms linked to 
them. The plasmonic projection (|37[) is now: 



d^ A + ^ 



whereas the soliton projection reads (|30|) 



-JJ2.A + P A 



+ N S 



dz 



;C + P 2 S C 



-A PP A 



A ps C 



-A K [2| A\ 2 C + A 2 C*] - T K [2\C\ 2 A + C 2 A*] 



Ag S C AgpA 



-A' K [2\A\ 2 C + A 2 C*} - T' K [2\C\ 2 A + C 2 A*} , 
where we have introduced two new Kerr coupling coefficients: 



r' 



K / 7 t np(0) fn P (z) j ? (Z) 



lfs(0)f ap (z)ff( Z ) 



However, since Ae p = 0, we know that their two projections vanish, so that, A sp 



and A z 



0. If we follow 



now the demonstration given in Section [V] we can conclude that the form of the generalized equation P^) has to be 
substituted by a new one including the new Kerr-coupling terms and in which q sp = (from its definition Eq. (|43[) 
along with A sp — and A pp = 0): 



dz 2 ' 

&_ 

dz 2 



Bl p )A = -qp S C-q K {2\A\ 2 C + A 2 C*)-pK(2\C\ 2 A + C 2 A* 
-q' K (2\A\ 2 C + A 2 C*) - p' K (2\C\ 2 A + C 2 A*) , 



ip 

B 2 I C 



(54) 



where we have introduced the new nonlinear coefficients 



Pk 



Pk 



n s t k 


^ps^K 


N ap N s 


— Sp S Ssp 




— Ssp^K 


N np N s 


SpsSsp 



(55) 
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We see that, as expected, despite there is no trace of linear plasmon-to-soliton coupling, nonlinear coupling occurs in 
the next-to-leading order weak coupling approximation. It is revealing to consider the case of a very weak plasmon 
field, in which we are approaching the linear plasmon limit and thus terms 0(A 2 ) can be neglected in Eq. (f54|) : 

q ps C~ PK (2\C\ 2 A + C 2 A*) 

P ' K (2\C\ 2 A + C 2 A*). (56) 

The previous equation makes clear that there is a second mechanism for plasmon-to-soliton coupling absent at leading 
order. Soliton can be also pumped or depleted by the plasmon field through a mechanism based on the Kerr coupling 
terms in Eq. (|56p . These terms can be understood as nonlinear sources generating gain or loss on the soliton param- 
eter C and driven by the plasmon field. As before, this mechanism can be visualized by renormalizing the soliton 
propagation constant in the following way: 



dz 2 ' 
dz 2 



A 



B 2 C 



B 2 =B 2 +p' K (2CTA + CA*), 
which permits to write the soliton equation as: 

d 2 



Analogously as in a previous analysis, the imaginary part of the renormalized propagation constant give us essential 
information about this nonlinear mechanism 

%(B 2 s )=p' K \C\\A\ sin (<p p -<p s ), 

which indicates that this type of term generates nonlinear gain or loss depending on the value of the soliplasmon 
relative phase. 

Another way of understanding the physical nature of the plasmon-to-soliton nonlinear term (the one proportional 
to p' K ) in Eq. (|56j) . and more specifically the one depending on |C| 2 , is by considering that its origin is the nonlinear 
modulation of the dielectric function in the Kerr medium. We know that the absence of linear modulation (Ae p = 0) 
makes the plasmon-to-soliton linear coupling to vanish. However, even if there is no linear modulation in the region 
where the soliton is localized (we assume an homogeneous nonlinear dielectric medium), the Kerr effect does induce 
a nonlinear modulation of the refractive index or, equivalently, of the dielectric function e{x) = Sk + 'y | i? s (a;) | 2 . 
Therefore, there exists, in fact, an effective modulation of s(x) induced by the Kerr nonlinearity that, in turn, 
originates a local variation of the dielectric function for the plasmon field, as defined in Eq. pTj) . given by 

Ae^ L {x) = e(x) - e p {x) = e(x) - e K = l\E s \ 2 if x > 0. 

Since this term is nonzero, an effective plasmon-to-soliton variational coupling, analogous to the linear one A sp 
(Ea.([4ip). is expected (recall that the linear coefficient vanish, A sp = 0) 

A^ = k 2 [ f s (0)Ae^fn P (z)^k 2 \C\ 2 [ 7f s (0)f!(z)f np (z)=T' K \C\ 2 . 
Jr. Jr 

The previous variational term corresponds to the soliton projection of Ae^ 1 ". One recognizes that this term is 
responsible of the first term of the p' K coefficient in Eq. (|55j) . This term can be interpreted as the leading order effect 
of soliton-soliton interaction in the weak coupling regime. In this case, the nonlinear plasmon, acting as a soliton, 
couples to the soliton tail though its own asymptotic exponential tail by means of the Kerr term. The soliton-soliton 
interaction exists independently of the existence of a linear modulation and, therefore, occurs even in a completely 
homogeneous medium. For this reason, it is expected to exist even if the soliton moves in a completely homogeneous 
medium, as in the present case in which we deal with a MK structure. Analogously, the nonlinear self-modulated 
local variation of the dielectric function Ae p L generates also a next-to-leading order coupling through its plasmon 
projection, which physically can be interpreted as generated by a nonlinearly induced term A^p of the type given in 
Eq.® (recall that there is no linear counterpart of this term since A pp — 0): 

. NL — 7,2 / ( (n s a ,-NL , 



A i ; p ^k 2 / f np (0)Ae^f np (z) = k 2 \C\ 2 / 7 /np(0)/ 2 (z)/„p(z)=r K |C| 



This coupling coefficient generates the second term in the expression for p' K (|39|) . 
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VII. CONCLUSIONS 

We have introduced a variational approach to properly understand the physics behind the problem of the nonlinear 
excitation of a SPP by a spatial soliton. Unlike in the original proposal in Ref.[l[, the resulting nonlinear resonator 
model is obtained from first principles, which permits, in this way, to provide an approximate solution of the full vector 
Maxwell's wave equation fl} for configurations close to our soliplasmon ansatz (f2"Tl) . In physical terms, variational 
equations extract the most relevant information of soliplasmon resonances as bound states of a soliton and a linear or 
nonlinear SPP. They provide the dynamics of a SPP and a soliton propagating along a metal/dielectric interface in 
the presence of a continuous exchange of electromagnetic energy between them in a process controlled by a nonlinear 
coupling constant, which is proportional to the soliton field at the interface e~ K " a (see Eq. (|50p ). Since the inverse 

1 /2 

soliton penetration length n s — (7/2) fco|C(z)| depends nonlinearly on the soliton amplitude, even the simpler case 
of a soliton bounded to a linear SPP presents special features with respect other well-know nonlinear resonator models 

On the other hand, the fact that soliplasmon resonances exists even in the presence of a low power SPP component, 
i.e, a linear SPP, demonstrates that nonlinearities in a soliplasmon play also a supplementary role than the one that 
is commonly attributed to them in nonlinear plasmonics. The fact that SPP's can support very high intensities very 
close to the metal/nonlinear-dielectric interface is the origin of most of the nonlinear effects reported in the literature. 
They are responsible of generation of second harmonic (see, e.g.. fl8l| and references therein) but also of nonlinearities 
of the Kerr type. They are known as early as in the 80's [3H1C|. More recently, substantial advances in the field 
of plasmonics have refreshed the possibility of exploiting these effects for nanophotonics applications using available 
technology, so a renovated interest in nonlinear effects induced by Kerr materials interacting with metals have been 
reflected in the literature (see, for example, (20l - [28| ). The high intensities reached in the dielectric in the vicinity 
of the metal interface associated to a high-power SPP mode are able to enhance the nonlinear response of a Kerr 
medium directly attached to the metal. This response can, in return, modify the propagation properties of the SPP. 
As pointed out as early as in Ref.[9], "when nonlinear dielectric media are in contact with a metal surface, the surface 
plasmon polaritons guided by that interface become power-dependent." This can be considered a standard definition 
of a nonlinear plasmon. The fact that the propagation constant of the SPP becomes power-dependent explicitly 
shows that, according to this definition of a nonlinear plasmon, the latter is continuously connected to the linear 
SPP in a standard P vs /i representation. From this perspective, soliplasmons are not nonlinear plasmons since their 
resonant nature forces their two branches — corresponding to 0- and n- soliplasmon solutions — to be detached from 
that of a linear SPP in a P vs /i representation (see its resonant behavior in Ref.[2f). Physically, in a soliplasmon, 
plasmon and soliton preserve their identity as localized solutions even though they are interacting. In a nonlinear 
plasmon, the soliton cannot be resolved as a second spatially localized component. Our variational approach for 
a soliplasmon formalizes this features explicitly by assuming the existence of two separately localized plasmon and 
soliton components in our ansatz ([2]) . In this way, our variational equations distinguish between two different types of 
nonlinear effects: (i) those affecting the uncoupled propagation of the plasmon and soliton components, i.e, those which 
appear even when they do not interact; and (ii) those related to the interaction. In the case of the SPP, the former 
are taken into account in our variational formalism through the dependence of the plasmon propagation constant 
Pnp = Pnp(\A\ 2 ) on the SPP amplitude. This is analogous to the well-known behavior of the soliton propagation 
constant with respect to its amplitude: f$\ = {3 2 (\C\ 2 ). The second type of nonlinear effects are related to coupling. 
Here we can, in turn, distinguish two different categories. The first one is related to typical soliton-to-soliton coupling, 
as those appearing in Eqs. (|51[) or (1561) . which show a form analogous to cross-phase modulation [19]. The second 
one is related to the modulation of the linear dielectric function and it is the analogous to the coupling between 
neighboring atoms in a crystal in the so-called tight binding approximation [14]. The coupling here is giving by the 
overlapping of the wave functions of two solutions of individual atomic sites detached from one another weighted by 
the difference of the total periodic potential with respect the atomic one. This is the nature of the coupling terms A sp 
and A ps involving the local variations of the linear dielectric function Ae p and Ae s in Eqs. ([38|) and (|41"j) . This is in 
fact the dominant coupling term in the weak coupling approximation, as demonstrated by our variational equations 
for the interaction of a linear plasmon and a soliton (|45[) . However, since the overlapping solutions are in this case 
nonlinear (in this case, that of the the soliton), the coupling becomes itself nonlinear. It is precisely this second type 
of nonlinearity in the coupling what causes the soliplasmon variational model to be qualitatively different from other 
coupled nonlinear systems |17j . 
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